/* 
 * File:   tersoffutils.hpp
 * Author: mik
 *
 * Created on July 13, 2013, 4:05 PM
 */

#ifndef TERSOFFUTILS_HPP
#define	TERSOFFUTILS_HPP
#include <math.h>
#include <assert.h>
#include <lattice/lattice.hpp>
#include <lattice/latticeutils.hpp>
#include <util/common/util.hpp>

#include "tersoffparams.hpp"

/**
 * Utils to compute Tersoff potential
 */
class TersoffUtils {
public:

    /**
     * Constructor
     * @param tpar pointer to the Tersoff parameters
     */
    TersoffUtils(TersoffParams* tpar) :
    mLatBase(NULL), mEps(0.0), mDelta(0.00001), mRange(2), mTpar(tpar) {
    }

    /**
     * Setup base lattice
     * @param lat pointer to the base (unscaled lattice)
     */
    void setLatticeBase(const Lattice* lat) {
        mLatBase = lat;
    }

    /**
     * Sets eps for numerical derivatives
     * @param eps to set
     */
    void setEps(double eps) {
        mEps = eps;
    }

    /**
     * Compute energy per atom
     * @return energy per atom
     */
    double energyPerAtom() {
        double v = computeEnergy();
        return v / mLatBase->mNCell;
    }

    /**
     * Compute cell energy 
     * @return computed energy
     */
    double computeEnergy() {
        return energy(mLatBase);
    }

    /**
     * Computes bulk module
     * @return bulk module computed
     */
    double computeBulk() {
        assert(mEps > 0);
        double e0 = energy(mLatBase) / (double) mLatBase->mNCell;
        Lattice tmp;
        LatticeUtils::copy(mLatBase, &tmp);

        double scl[3];
        double a = 1 + mEps;
        scl[0] = a;
        scl[1] = a;
        scl[2] = a;
        LatticeUtils::scale(scl, mLatBase, &tmp);
        double e1 = energy(&tmp) / (double) mLatBase->mNCell;

        a = 1 - mEps;
        scl[0] = a;
        scl[1] = a;
        scl[2] = a;
        LatticeUtils::scale(scl, mLatBase, &tmp);
        double e2 = energy(&tmp) / (double) mLatBase->mNCell;

        double c = 1.602 * 8. / (9. * mLatBase->mLC);
        double u = c * (e1 + e2 - 2 * e0);
        double dx = mLatBase->mLC * mEps;
        double B = u / (dx * dx);

        return B;
    }

    /**
     * Computes shear module
     * @return shear module computed
     */
    double computeShear() {
        assert(mEps > 0);
        double e0 = energy(mLatBase) / (double) mLatBase->mNCell;
        Lattice tmp;
        LatticeUtils::copy(mLatBase, &tmp);

        double scl[3];
        scl[0] = 1 + mEps;
        scl[1] = 1. / (1. + mEps);
        scl[2] = 1;
        LatticeUtils::scale(scl, mLatBase, &tmp);
        double e1 = energy(&tmp) / (double) mLatBase->mNCell;

        scl[0] = 1 - mEps;
        scl[1] = 1. / (1. - mEps);
        scl[2] = 1;
        LatticeUtils::scale(scl, mLatBase, &tmp);
        double e2 = energy(&tmp) / (double) mLatBase->mNCell;

        double c = 1.602 * 2. / mLatBase->mLC;
        double u = c * (e1 + e2 - 2 * e0);
        double dx = mLatBase->mLC * mEps;
        double C = u / (dx * dx);

        return C;
    }

    /**
     * Computes C4 module and Kliemann constant
     * @param kli
     * @return 
     */
    double computeC44Base(double& kli) {
        assert(mEps > 0);
        double e0 = energy(mLatBase) / (double) mLatBase->mNCell;
        Lattice tmp, tt;
        LatticeUtils::copy(mLatBase, &tmp);
        LatticeUtils::copy(mLatBase, &tt);
        double dx = 0.00001;

        double tr1[3][3] = {
            {1, 0.0, 0},
            {0, 1, 0},
            {0, 0, 1}
        };
        tr1[0][1] = dx;
        LatticeUtils::transform((double*) tr1, mLatBase, &tmp);
        double e1 = energy(&tmp) / (double) mLatBase->mNCell;
        double d = 0.0001;
        double vec[3] = {0, 0, 0};
        double emin = e1;
#define TERSUT_COMPUTE_OBJECTIVE(X, E)\
{\
vec[2] = X;\
LatticeUtils::shiftInternal(vec, &tmp, &tt);\
E = energy(&tt) / (double) mLatBase->mNCell;\
}
        double xl = -d, xr = d;
        double fl, fr;
        TERSUT_COMPUTE_OBJECTIVE(xl, fl)
        TERSUT_COMPUTE_OBJECTIVE(xr, fr)
#if 1        
                int n = 0;
        while (xr - xl >= 0.000001) {
            double l = xr - xl;
            double x1 = xl + l / 3;
            double x2 = xr - l / 3;
            double f1, f2;
            TERSUT_COMPUTE_OBJECTIVE(x1, f1)
            TERSUT_COMPUTE_OBJECTIVE(x2, f2)
            if (f1 > fl) {
                xr = x1;
                fr = f1;
            } else if (f2 > fr) {
                xl = x2;
                fl = f2;
            } else if (f1 > f2) {
                xl = x1;
                fl = f1;
            } else {
                xr = x2;
                fr = f2;
            }
            n++;
        }
        e1 = 0.5 * (fl + fr);
        double dz = BNBABS(0.5 * (xl + xr));
        kli = 4. * dz / (dx * mLatBase->mLC);
#endif 
#if 0
        printf("emin = %.24lf, \ne1 = %.24lf\n", emin, e1);
        e1 = emin;

        double tr2[3][3] = {
            {1, -0.00001, 0},
            {0, 1, 0},
            {0, 0, 1}
        };
        LatticeUtils::transform((double*) tr2, mLatBase, &tmp);
        double e2 = energy(&tmp) / (double) mLatBase->mNCell;
        printf("e2 = %.24lf\n", e2);
#endif
        double e2 = e1;

        double c = 1.602 * 8. / pow(mLatBase->mLC, 3);
        double u = c * (e1 + e2 - 2 * e0);
        double C = u / (dx * dx);

        return C;
    }

    /**
     * Computes lattice constant as the argmin of potential energy
     * @return lattice constant
     */
    double computeLatticeConst() {
        assert(mDelta > 0);
        double e0 = energy(mLatBase);
        double minlc;
        double base = mLatBase->mLC;
        double xl = base - mRange, xr = base + mRange;
        double fl, fr;
        fl = scaledEnergy(xl / base);
        fr = scaledEnergy(xr / base);
        int n = 0;
        while (xr - xl >= mDelta) {
            double l = xr - xl;
            double x1 = xl + l / 3;
            double f1 = scaledEnergy(x1 / base);
            double x2 = xr - l / 3;
            double f2 = scaledEnergy(x2 / base);
            if (f1 > fl) {
                xr = x1;
                fr = f1;
            } else if (f2 > fr) {
                xl = x2;
                fl = f2;
            } else if (f1 > f2) {
                xl = x1;
                fl = f1;
            } else {
                xr = x2;
                fr = f2;
            }
            n++;
        }

        double rv = 0.5 * (xl + xr);
        return rv;
    }

    void printValues() {
        printf("Energy per atom: %lf\n", energyPerAtom());
        printf("Bulk module: %lf\n", computeBulk());
        printf("Compute shear: %lf\n", computeShear());
        double k;
        double c44 = computeC44Base(k);
        printf("C44 module: %lf\n", c44);
        printf("Kliemann: %lf\n", k);
        printf("Lattice const: %lf\n", computeLatticeConst());
    }
    
    /**
     * Retrieve tersoff parameters
     * @return Tersoff parameters
     */
    TersoffParams* getTpar() {
        return mTpar;
    }
private:

    double scaledEnergy(double s) {
        Lattice tmp;
        LatticeUtils::copy(mLatBase, &tmp);
        double scl[3];
        scl[0] = s;
        scl[1] = s;
        scl[2] = s;
        LatticeUtils::scale(scl, mLatBase, &tmp);
        double e = energy(&tmp);
        return e;
    }

    /**
     * Compute cell energy 
     * @return computed energy
     */
    double energy(const Lattice* lat) {
        assert(lat);
        mLat = lat;
        double rv = 0;
        for (int i = 0; i < mLat->mNCell; i++) {
            double ei = 0;
            for (int j = 1; j < mLat->mNei[i][0]; j++) {
                double eij = computePairEnergy(i, mLat->mNei[i][j]);
                ei += eij;
            }
            rv += ei;
        }
        return 0.5 * rv;
    }

    /**
     * Cutoff function
     * @param r radius
     * @return cut value between 0 and 1
     */
    double cutoff(double r) {
        double rv;
        if (r < mTpar->mR - mTpar->mRcut) {
            rv = 1;
        } else if (r > mTpar->mR + mTpar->mRcut) {
            rv = 0;
        } else {
            rv = 0.5 * (1 - sin(M_PI * (r - mTpar->mR) / (2. * mTpar->mRcut)));

        }
        return rv;
    }

    /**
     * Repulsive energy
     * @param r distance
     * @return repulsive energy
     */
    double VR(double r) {
        double rv;
        rv = mTpar->mDe / (mTpar->mS - 1.);
        rv *= exp(-mTpar->mB * sqrt(2. * mTpar->mS) * (r - mTpar->mRe));
        return rv;
    }

    /**
     * Attractive energy
     * @param r distance
     * @return attractive energy
     */
    double VA(double r) {
        double rv;
        rv = mTpar->mDe * mTpar->mS / (mTpar->mS - 1.);
        rv *= exp(-mTpar->mB * sqrt(2. / mTpar->mS) * (r - mTpar->mRe));
        return rv;
    }

    /**
     * Compute energy of two atoms interaction
     * @param i the number of the first atom
     * @param j the number of the second atom
     * @return the computed energy
     */
    double computePairEnergy(int i, int j) {
        double r = computeDist(i, j);
        double rv = cutoff(r);
        if (rv > 0) {
            rv *= (VR(r) - computeB(i, j) * VA(r));
        }
        return rv;
    }

    double computeDist(int i, int j) {
        double rv = 0;
        for (int k = 0; k < 3; k++)
            rv += (mLat->mCoor[i][k] - mLat->mCoor[j][k])*(mLat->mCoor[i][k] - mLat->mCoor[j][k]);
        return sqrt(rv);
    }

    double computeB(int i, int j) {
        double rv;
        double zeta = computeZeta(i, j);
        double a = 1. + pow(mTpar->mGamma * zeta, mTpar->mNu);
        rv = pow(a, -1. / (2. * mTpar->mNu));
        return rv;
    }

    double computeZeta(int i, int j) {
        double rv = 0;
        double rij = computeDist(i, j);
        for (int k = 0; k < mLat->mNA; k++) {
            if ((k != i) && (k != j)) {
                double rik = computeDist(i, k);
                double cut = cutoff(rik);
                if (cut > 0) {
                    double e;
                    double v = mTpar->mLambda * (rij - rik);
                    double w = exp(v * v * v);
                    e = cut * w * computeG(i, j, k);
                    rv += e;
                }
            }
        }
        return rv;
    }

    double computeG(int i, int j, int k) {
        double rv = 1;
        double v = mTpar->mC / mTpar->mD;
        rv += v * v;
        double *xi = (double*) (mLat->mCoor[i]);
        double *xj = (double*) (mLat->mCoor[j]);
        double *xk = (double*) (mLat->mCoor[k]);
        double s = (xj[0] - xi[0])*(xk[0] - xi[0]) + (xj[1] - xi[1])*(xk[1] - xi[1]) + (xj[2] - xi[2])*(xk[2] - xi[2]);
        double rij = computeDist(i, j);
        double rik = computeDist(i, k);
        double c = s / (rij * rik);
        v = mTpar->mH - c;
        rv -= mTpar->mC * mTpar->mC / (mTpar->mD * mTpar->mD + v * v);
        return rv;
    }


    TersoffParams* mTpar;
    const Lattice* mLatBase;
    const Lattice* mLat;
    double mEps;
    double mDelta;
    double mRange;
};

#endif	/* TERSOFFUTILS_HPP */

